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This work is concerned with the Biot/squirt (BISQ) models for wave prop¬ 
agation in saturated porous media. We show that the models allow exponen¬ 
tially exploding solutions, as time goes to infinity, when the characteristic 
squirt-flow coefficient is negative or has a non-zero imaginary part. We also 
show that the squirt-flow coefficient does have non-zero imaginary parts for 
some experimental parameters. Because the models are linear, the existence 
of such exploding solutions indicates instability of the BISQ models. This 
result calls on a reconsideration of the widely used BISQ theory. Further¬ 
more, we demonstrate that the 3D isotropic BISQ model is stable when the 
squirt-flow coefficient is positive. In particular, the original Biot model is un¬ 
conditionally stable where the squirt-flow coefficient is 1. 
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1. Introduction 

The propagation of seismic waves in rocks with fluids has been among the most active 
research fields in geoscience. Various theories [Biot, 1956a, b; Dvorkin et ah, 1993; Mavko 
et ah, 2009; Muller et al, 2010] focus on relating attenuation and dispersion of seismic 
waves in earth materials to physical properties of the rocks and fluids. Attenuation means 
the exponential decay of wave amplitude with distance and dispersion is a variation of 
propagation velocity with frequency. Among the various mechanisms related to attenu¬ 
ation and dispersion, the Biot and squirt-flow mechanisms are believed to be the most 
important ones [Dvorkin et al., 1993; Yang et al., 2002], They have served as rigorous 
and formal foundations to study acoustic wave propagation in saturated porous media. 

In Biot’s theory [Biot, 1956a, b], pore fluids are forced to participate in the solid’s motion 
due to inertia and viscosity. The Biot mechanism is described through macroscopic rock 
properties and based on the macroscopic flow. The squirt-flow mechanism is associated 
with the squirting of the pore fluid out of cracks as they are deformed by the compressional 
waves. It is related to microscopic rock properties and based on the local fluid flow 
[Dvorkin et al, 1993]. Clearly, the two mechanisms occur simultaneously and affect the 
process of seismic energy propagation and attenuation. 

Dvorkin and Nur [1993] developed a consistent one-dimensional (ID) model, incorporat¬ 
ing both the Biot and squirt-flow mechanisms. That is the widely used Biot/Squirt-flow 
(BISQ) model. The BISQ theory directly relates the velocity and attenuation of P-waves 
to the macroscopic and microscopic parameters. Later on, Parra [1997, 2000] extended 
the ID model to the 3D case for a transversally isotropic poroelastic medium, which is 
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based on the constitutive relations in the frequency domain given by Kazi-Aoual et al. 
[1988] and the isotropic description of the solid/fluid coupling effects developed by Biot 
[1956b], Furthermore, Yang and Zhang [2002] adapted the BISQ theory for anisotropic 
porous media. Up to now, it is well accepted [Mavko et al, 2009; Midler et al, 2010] that 
the BISQ theory is adequate to model acoustic wave propagation in saturated porous me¬ 
dia. The acceptance is based mainly on experimental verifications of the predicted phase 
velocities or analysis of wave field snapshots for a short time. However, a mathematical 
study of the BISQ model does not seem available in the literature. 

Observe that the models can be transformed to a system of first-order partial differ¬ 
ential equations (see Section 3) and thereby possess the form of hyperbolic relaxation 
systems. Based on this, we are guided by the theory of hyperbolic relaxation systems 
developed in [Yong, 1992, 2001] and conduct a stability analysis of the BISQ model and 
its 3D isotropic extension. Our results show that the models allow spatially bounded but 
time exponentially exploding solutions when the characteristic squirt-flow coefficient is 
negative or has a non-zero imaginary part. By numerically examining some experimental 
parameters from [Dvorkin et al, 1993; King et ah, 2000; Parra, 2000; Yang et ah, 2011], 
we see that it is usual and realistic for the squirt-flow coefficient to have non-zero imag¬ 
inary parts and for the models to have time exponentially exploding solutions. Because 
the models are linear, the existence of such exploding solutions indicates instability of 
the widely used models. This seems the first time to point out instability of the BISQ 
models. Furthermore, we show that, when the squirt-flow coefficient is positive, the 3D 
isotropic model satisfies a structural stability condition for hyperbolic relaxation systems 


DRAFT 


May 12, 2015, 3:19am 


DRAFT 



LIU ET AL.: STABILITY ANALYSIS OF THE BIOT/SQUIRT MODELS 


X - 5 


and therefore is stable. In particular, the original Biot model is unconditionally stable 
where the squirt-flow coefficient is 1. 

Our results call on a reconsideration of the BISQ theory. A possible modification 
would be to replace the squirt-flow coefficient with its real part. The results may also 
provide reasonable explanations for some problematic numerical simulations reported in 
the literature. In [Parra, 2000] for a transversely isotropic medium, it was pointed out that 
the velocity dispersion and attenuation predicted by the BISQ model are not consistent 
with experimental data at the low angular frequency. For the snapshot of seismic wave 
fields, Yang et al. [2002] found that it is hard to gain the slow P-waves by the BISQ model 
with the small frequency and the slow P-waves occur only after increasing the frequency. 
In addition, the numerical results in Yang et al. [2002] indicate that the slow P-waves only 
occur when the fluid viscosity is negligible. On the other hand, we notice that the squirt- 
flow coefficient has always a nonzero imaginary part at small angular frequencies and is 
real when the viscosity vanishes. These coincide very well with our theoretical results: 
the BISQ model is unstable when the squirt-flow coefficient has a nonzero imaginary part 
and is stable when the squirt-flow coefficient is positive. 

The paper is organized as follows. In Section 2 we show the existence of tirne- 
exponentially exploding solutions to the ID BISQ model when the squirt-flow coefficient 
is negative or has a non-zero imaginary part. Here we also numerically examine some 
experimental parameters. Section 3 is devoted to the 3D isotropic extension, which is 
transformed to a system of first-order partial differential equations with the coefficient 
matrices given in Appendix. It is demonstrated that the instability result is valid for the 
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3D model and the positive squirt-flow coefficient ensures the structural stability of the 
extension. The final section summarizes our main results. 

2. ID Model 

We start with the Biot/squirt (BISQ) dynamic equations for a two-component solid- 
fluid continuum in one dimension [Dvorkin et ah, 1993]: 

(1 - (ft)p s u a + < i>pfW tt = Mu xx - aP x , 

4>pfiv a - p a (u tt - w tt ) - [pcft 2 (u t - W t )]/K = ~<j)P x , (1) 

Pt = ~FS {w xt + [(a - (j))u xt }/(j)} . 

Here u = u(x,t ) and w = w(x,t ) are the respective solid’s and fluid’s displacements at 
the position-time point (x, f), P = P(x,t ) is the total fluid pressure, (ft G (0,1) is the 
porosity of the solid, M is the uniaxial modulus of the skeleton in drained conditions, a 
is the poroelastic coefficient, p s , pj and p a are the respective solid’s, fluid’s and additional 
coupling densities, p is the viscosity of the fluid, k is the permeability of the skeleton, 
and F is the Biot-flow coefficient. All these parameters are positive (see Table 1 below). 
Subscripts x and t indicate partial derivatives ( P x — -|j, u tt — etc.). 

In (1), S is the characteristic squirt-flow coefficient, ft has the expression 

5 = 1 — [2 Ji (AR) ] / [Ai? J 0 (A-R) ] (2) 

with 

A2 = {Pff 2 [ l + Pa/{<t>Pf ) + V-li:r](ft)/{np f f)] } /F, (3) 
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where Jo and J\ the respective Bessel functions of zero- and first-order, R is the charac¬ 
teristic squirt-flow length, and / is the angular frequency. Observe that A and thereby S 
are complex in general. See also Table 1 below. 

By substituting the third equation in (1) into the first two equations therein, we get 
two equations for the displacements u and w. For the resultant two equations, we look 
for their solutions of the form [Yong, 1992] 

u(x,t ) = u(t)exp(i£x),w(x,t) = w(t)exp(i£x) 
with £ real. In this way, we obtain 

PiUtt + p- 2 Wu = - {M + [oiFS(a - 0)]/0}£ 2 w - aFS^w, 

pi 2 U tt + P 22 W tt ~ [' ri4> 2 {u t - w t )\/n = (f)FS£ 2 {-[(a - (j>)u\/<j> - w) . (4) 

Here pi = (1 - <j>)p a , p 2 = <j>p f , pvi = ~p a and p 22 = </>p f + p a - 

Set v = u t and V = w t . The system of ordinary differential equations (3) can be 
rewritten as 

u t = v, 

w t = V, 

Piv t + p 2 Vt = - {M + [aFS(a - 4>)\/4>} £ 2 n - aFS£ 2 w, 

pi 2 V t + P 22 Vt = -i 2 FS(a - 4>)u - £ 2 (f)FSw + (r)(f) 2 v)/K - (p0 2 H)/k. 
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This can be further rewritten as 


AU t = BU. 


( 5 ) 


Here U = [u,w,v,V] T is an unknown vector-valued function of £ € [0,+oo), A and B 


are constant coefficient matrices 
10 0 0 
0 10 0 
0 0 pi p 2 
0 0 pi 2 p 2 2 


A = 


B = 


0 0 1 

0 0 0 

-£ 2 [M + aFS(a- <t>)/<j>)} -£, 2 aFS 0 


0 

1 

0 


-eFS(a - 0) 


-^(pFS (i/0 2 )/k -(p0 2 )/k 


The solution to the system (5) with initial data Uq is 


U(t) = exp(tA _1 B)U 0 . 


This solution will exponentially explode as t goes to infinity, provided that the matrix 
A V B has an eigenvalue with positive real part and Uo is the corresponding eigenvector. 
Because the BISQ model (1) is linear, the existence of such exploding solutions means 
instability. 

The eigenvalues of matrix A l B are zeros of its characteristic polynomial 

A 4 + (p&A 3 ) / k + £ 2 cA 2 + (p^ 2 dA) / k + £ 4 e = 0 (6) 


with coefficients 


b = [0 2 (pi + p 2 )]/(pip 2 2 - P 2 P 12 ), 


C = {M0P22 + [(pi + p 2 )0 2 - 2p 2 «0 + p220l 2 ]FS) / [0(piP22 ~ p2pn)\ , 

d=y> (M0 + a 2 F5)] /(P 1 P 22 - P 2 P 12 ), e = (M0FS')/(pip 22 - p 2 pu)- (7) 
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These coefficients are derived through a direct calculation and also with Mathlab. Notice 
that (P 1 P 22 — P2P12) is always positive for p 12 < 0. 

About the matrix A l B, we have 

Proposition 2.1 If S < 0, then the matrix A l B withf ^ 0 has eigenvalues with positive 
real parts. 

Proof. In case S is real, the coefficients of the characteristic polynomial are all real. 
Thus, the complex eigenvalues must be conjugate. If S < 0, the product of four eigenvalues 
£ 4 e is negative. This means that there exists at least two real eigenvalues with opposite 
sign. Hence the matrix A ~ l B has eigenvalues with positive real parts. 

Recall that the squirt-flow coefficient S is complex in general. We point out the following 
important fact. In what follows, we denote by Re(a) and Im(a) the real and imaginary 
parts of complex number a, respectively. 

Theorem 2.1 If the imaginary part of S is non-zero, then the matrix A L B has eigen¬ 
values with positive real parts whenever fn/g is sufficiently large (positive or negative). 
Proof. Observe that the roots of the algebraic equation (6) can be written as 

A = £A [^/(f/c)] 

with A = A(e) solving 

A 4 + eAA 3 + cA 2 + edX + e = 0. 

Note that A(e) is continuous with respect to e [Kato, 1982], Then p = lim^o A(e) solves 

p 4 + cp 2 + e = 0. 


DRAFT 


May 12, 2015, 3:19am 


DRAFT 



X - 10 LIU ET AL.: STABILITY ANALYSIS OF THE BIOT/SQUIRT MODELS 

Denote by p 3 (j = 1,2, 3,4) the roots of this algebraic equation. Because their product e 
is complex whenever Im(S)^ 0, there must be a j\ e {1,2, 3,4} such that Re(/ij 1 ) 0. 

Since Jfj Pj = there must be another j 2 G {1,2, 3,4} such that Re(/ij 1 )Re(/ij 2 ) < 0. 
Thus, by the continuity argument the algebraic equation (6) has eigenvalues with positive 
real parts for large /r). This completes the proof. 

This theorem can be strengthened as follows. 

Theorem 2.2 Under the condition 



the matrix A l B has eigenvalues with positive real parts. 

Proof. Let A j = Xj + \/—ly 3 (j = 1,2, 3,4) be solutions to the algebraic equation (6). 

Notice that b = [ 4 > 2 (pi + p2)\/(piP22 — P2P12) > 0. We have Yh x j = —1 ^ >r l/ K i = 0 and 

j j 

thereby -2 Y, ]<k VjVk = yl + vl + vl + hi 
If Xj < 0 for all j, we deduce that 

2£ Re(c) 2 Re ^ ^ AjA^ 2 XjX k DjUk ) 2 ^ ^ x 3 x k + ^ ^ yj U ^ ^ Vj 

j<k j<k j<k j j 

and 
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The above results indicate that the ID BISQ model (1) allows spatially bounded but 
time-exponentially exploding solutions when the squirt-flow coefficient is negative or has 
a non-zero imaginary part. Theorem 2.2 gives a sufficient condition for the existence of 
the exploding solutions. 

Motivated by the above results, we compute the eigenvalues with parameters from Refs. 
[Dvorkin et al, 1993; King et al, 2000; Parra, 2000; Yang et al., 2011] and £ = 1. The 
following table contains the imaginary parts of FS and the corresponding eigenvalues 
with the largest real parts, which are all positive. 


Table 1. Eigenvalues with experimental parameters and £ = 1 


Data Sources 

[Dvorkin et al., 1993] 

[King et al., 2000] 

[Parra, 2000] 

[Yang et al., 2011] 

M{GPa) 

35.480 

20.217 

23.278 

29 

0 (%) 

15 

27.4 

23 

15 

a 

0.5789 

0.7692 

0.6197 

0.5702 

p s (kg/m 3 ) 

2650 

2550 

2750 

2600 

p f (kg/m 3 ) 

1000 

1000 

1000 

1000 

p a (kg/m 3 ) 

420 

420 

420 

450 

77 /k ( kg/m 4 • s) 

4 x 10 4 

5xl0 4 

8 xl 0 4 

100 

/ (KHz) 

50 

600 

50 

50 

R (m) 

0.001 

0.015 

0.001 

0.001 

Im(FS) 

-37.5587 

-2.049X10 4 

-115.1 

-0.0939 

Max[Re(A)] 

16.9131 

6.018X10" 13 

14.7252 

17.6842 


3. 3D isotropic models 

The poroclastic wave equations for a 3-D motion including the Biot and squirt-flow 
mechanisms in an isotropic medium read as [Parra, 1997; Yang et al, 2002] 

c)~ 'ti) 

/iVV ■ u + pAu + {A + [aFS(a - 4>)\/<t>} VV • u + aFS VV • w = pi—r- + 

ot z ot z 

FS(a> - 0)VV • rr + (pFSW ■ w = p V2 ^ + p 22 + [# 2 - ^J]/«. (8) 
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Here A and /i are the Lame coefficients, u = u(xi,x 2 ,x 3 ,t) and w = w(xi,x 2 ,x 3 ,t) 
denote the respective solid’s and fluid’s displacement vectors, and the other parameters 
are same as those in (1) with M = A + 2 ft. 

Let [u(x,t),w(x,t)] solve the equations in (1). It is direct to verify that 

'ti rr 

u(x 1 ,x 2 ,x 3 ,t) = 0 u(xi,t) w(x 1 ,x 2 ,x 3 ,t) = 0 w(xi,t) 

0 J [o 

is a solution to equations (7). Thus, we have 

Proposition 3.1 If S < 0 or Im(S ) ^ 0, then equations (7) have spatially bounded but 
time-exponentially exploding solutions. 

Next we study the case S > 0, including the Biot model where S — 1. First of all, we 
transform the second-order equations (7) into a system of first-order partial differential 


equations. To do this, we write 



and set 
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for i , j = 1,2,3. With such notations, the linear system (7) can be rewritten as 
Pifri + ftSL - {A + [ a CS( a - 0)]/0} £(£ Wj) - 2p&W, - T n - u£-T a - a£L = 0, 

J = 1 

P.ft-2 + P2 - {A + [aFS(a - 0)]/0} gfj(E Wj) - 2 p£w 2 - - p£-T 23 - a£;L = 0, 

3 =1 

P.|« 3 + p 2 f U - {A + [ a cs( a - 0)]/0} ^(E - 2p^w 3 - p^r 13 - p£-t 23 - a £-p = o, 

3 = 1 

+ P2 2 |V - WS(a - </>)V(£ W,-) - 0VL = -[# 2 (V - «)]/«. (10) 

i=i 


In addition, by definition it is easy to see that 


_ UlE _ n Oljj _ Ovi_ _ SOU — 0(3 ^ 3\ OL pQ 

at, axj u > at, a.,- . a.r, U U ^ J)i at, 10 


din _ 


9T U 


dvj 


dvj 


dL 


dxi 


dxi 



= 0. 


( 11 ) 


In (9) and (11) there are 13 first-order partial differential equations. Clearly, they can be 
written in the vector form: 


4 u+ ^ A ’£-j u=Bu (i2) 

j =1 j 

with U = [W\, W 2 , W 3 , T 12 , T 13j T 23 , L, vi, Vi, f 2 , V 2 , V 3 , Vs] T and the coefficient matrices 
given in Appendix. Note that A is invertible. 

This system of first-order partial differential equations possesses the form of hyperbolic 
relaxation systems, whose fundamental properties were systematically identified in [Yong, 
2001]. Here we show that the system ( 12 ) with S > 0 satisfies the second structural 
stability condition in [Yong, 1992, 2001]. More precisely, we have 

Theorem 3.1 If S > 0, then there exists a symmetric positive definite matrix A 0 such 
that the matrices A T A 0 Aj (j = 1,2,3) are symmetric and A 1 A 0 B is symmetric non- 
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positive definite. Moreover, the L 2 -norm of the solutions to the system (12) are uniformly 

bounded in time. 

Proof. The symmetric positive definite matrix A 0 is given in Appendix. One can 
directly verify that A 1 Aq A j (j = 1,2,3) are symmetric and A r A$B is symmetric non¬ 
positive definite. 

To show the uniform boundedness, we multiply the system (12) with U 1 A 1 Aq and use 
the symmetry of A 7 AqAj (j = 1, 2, 3) to obtain 

4[ U t A t A 0 AU} + V -^—[U T A T A 0 AjU] = 2 U T A T A 0 BU. 

at ^ ax; 

j J 

Integrating this equation over [x\,x 2 ,xf\ G B? gives 

[ U 1 A 1 A Q AUdxidx 2 dx 3 = 2 / U T A T A 0 BUdxidx 2 dx 3 < 0, 

9t J R 3 J R 3 

where the inequality is clue to the nonpositive definiteness of the symmetric matrix 
A 1 A 0 B. Let the largest and smallest eigenvalues of the symmetric positive definite 
matrix A 7 AqA be A and a. Then it follows from the last inequality that 


Jr.U t (xi,x 2 , x 3 , t)U{x\, x 2 , x 3 , t)dx\dx 2 dx 3 

< a " 5& uT (x\,x 2 , x 3 , t)A r A 0 AU(x i, x 2 , x 3 , t)dx\dx 2 dx 3 

< a " !r> uT (xi,x 2 , x 3 , 0)A 3 A 0 AU(x i, x 2 , x 3 , 0)dxidx 2 dx 3 

< a~ l A Jr 3 U T (x 1 , x 2 , x 3 , 0)U(xi , x 2 , x 3 , 0)dxidx 2 dx3, 


namely, 


\\U(-,t)\\ 2 L 2 < a _1 A||?7(-, 0)||^2 


for all t > 0. This completes the proof. 
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The last theorem indicates the stability of the 3D isotropic BISQ model (7) with any 
parameters satisfying S > 0. Particularly, the original Biot model [Biot, 1956a, b] as the 
high-frequency limit of the 3D BISQ model (7) is unconditionally stable. On the other 
hand, we can directly verify the stability of the Biot model when R{2N + A) — Q 2 > 0 
(Biot’s parameters). The latter was identified also by Biot in [Biot, 1956a] with physical 
considerations, that is, it is equivalent to the positiveness of the kinetic energy. 

4. Conclusions 

In this work, we conduct a critical assessment of the BISQ models for isotropic poroe- 
lastic media and obtain the following conclusions. 

• When the squirt-flow coefficient S is positive, the 3D isotropic BISQ model is stable. 
Particularly, the original Biot model [Biot, 1956a, b] is unconditionally stable. 

• Otherwise, the model allows spatially bounded but time-exponentially exploding so¬ 
lutions and therefore are unstable. 

• Theorem 2.2 gives a sufficient condition for the existence of the exploding solutions, 
which requires that £k/t) is large enough when Im(S') ^ 0. Table 1 further shows that it 
is realistic and usual to have Im(S') ^ 0 and eigenvalues with positive real parts. 

These, for the first time, reveal the instability of the BISQ models and call on a recon¬ 
sideration of the widely used BISQ theory. 
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5. Appendix. Coefficient matrices for the 3D isotropic model 


In Equation (12), the coefficient matrices are 


B = -[(# 2 )/k 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 -1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-1 

1 J 


1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

o 1 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Pi 

P2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

P12 

P22 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Pi 

P2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

P12 

P22 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

pi 

P2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

P12 

P22 - 



r 0 

0 

0 

0 

0 

0 

0 

-1 


0 

0 0 

0 

0 1 


0 

0 

0 

0 

0 

0 

0 

0 


0 

0 0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0 

0 0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0 

-1 0 

0 

0 
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where A\ = 2/x + A 2 , A 2 = A + [aFS(a — (f>)]/(f>, A 3 = FS(a — 4>). 
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The matrix Aq in Theorem 3.1 is constructed as 

r A 4 A 5 A 5 0 0 0 Aq 0 0 0 0 0 0 

A 5 A 4 A5 0 0 0 Aq 0 0 0 0 0 0 

A 5 A5 A 4 0 0 0 Aq 0 0 0 0 0 0 

0 0 0 p/(pi+p 2 ) 00 0000000 

0 0 0 0 p/(pi+p 2 ) 0 0000000 

0 0 0 0 0 At/(pi+P 2 ) 0 0 0 0 0 0 0 

Ao = A e Aq Ae 0 0 0 <f>/[FS(pi+p 2 )] 0 0 0 0 0 0 

0000 0 0 0 As —A 7 0 0 0 0 

0 0 0 0 0 0 0 —A 7 (P 1 P 22 — P 2 P 12 )—* 0 0 0 0 

0000 0 0 0 0 0 As -At 0 0 

0000 0 0 0 0 0 —A 7 (pip 22 —P 2 P 12 ) ^ 0 0 

0000 0 0 0 0 0 0 0 As -A 7 

.0000 0 0 0 0 0 0 0 —At (pip 22 —p 2 pi 2) —1 

with 


A 4 — {A + 2 /i + [FS(a — 0) 2 ]/0}/(pi + P 2 ), A 5 — {A + [FS(a — <f>Y)/<t >}/(pi + P 2 ), 


Aq — (a — 4>)/(pi + P 2 ), A 7 — (p 12 + P 22 ) / [(pi + P 2 XP 1 P 22 — P 2 P 12 )], 


Ag — (pi + P 2 ) 2 + (P 12 + P 22 ) 2 / [(Pi + P 2) 2 (PlP 22 — P 2 P 12 )] • 


Recall that S > 0 and P 1 P 22 — P 2 P 12 > 0 for p 12 < 0. It can be directly checked that the 


matrix A 0 is symmetric and positive definite. 
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